#include "kernel/driver.hpp"

#include "yocto/string/vfs-utils.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/ios/ocstream.hpp"
#include <iomanip>

static const char fmt[]=" %.15f";

static inline 
void save_data( ios::ocstream &fp, Real t, const Cell &cell, const Real PCO2 )
{
	const Solution &S_in          = cell.solutions[ Cell::iS_in ];
	fp("%.4g",t/60.0);
	fp(fmt,1000.0*cell.Em);
	fp(" %.6f", S_in.pH());
	S_in.write_values(fp, fmt);
	fp(" %.6f", PCO2*760.0);
	fp("\n");
	
}


#if defined(SAVE_LEAK)
static inline
void save_leak( ios::ocstream &fp, Real t, Cell &cell )
{
	const Solution &S_in          = cell.solutions[ Cell::iS_in ];
	const Solution &S_out         = cell.solutions[ Cell::iS_out ];
	Solution       &dS            = cell.solutions[ Cell::idS    ];
	const Real      zeta_m        = cell.Em * cell.F_over_RT;
	fp("%.4g",t/60.0);
	fp(fmt,1000.0*cell.Em);
	fp(" %.6f", S_in.pH());
	dS.zero();
	dS.addLeak(zeta_m, cell.Volume, S_in, S_out);
	dS.write_values(fp, fmt);
	fp("\n");
}
#endif

int main( int argc, char *argv[] )
{
	
	const char *progname = _vfs::get_base_name( argv[0] );
	try
	{
		//======================================================================
		//
		// preparing data
		//
		//======================================================================
		if( argc <= 4 )
			throw exception("Usage: %s config.lua dt tmax output_name", progname);
		const string config_name = argv[1];
		
		
		Lua::State   VM;
		lua_State   *L = VM();
		lua_register(L,"geom_data",geom_data);
		Lua::Config::DoFile( L, config_name );
		Cell         cell( L );
		
		cell.chemsys.ftol = 1e-8;  // frational tolerance
		cell.chemsys.dt   = 0.1;   // P_CO2 time scale if needed
		
		Lua::Config::DoString(L, vformat("dt=%s;tmax=%s;", argv[2], argv[3]));
		const Real dt   = Lua::Config::Get<lua_Number>(L,"dt");
		const Real tmax = Lua::Config::Get<lua_Number>(L,"tmax");
		
		if( dt <= 0 )
			throw exception("Invalid dt=%g", dt );
		
		const string output_name = argv[4];
		
		//======================================================================
		//
		// initialize steady state
		//
		//======================================================================
		std::cerr << std::scientific;
		std::cerr.precision(15);
		const Solution &S_in          = cell.solutions[ Cell::iS_in ];
		Solution       &S_tmp         = cell.solutions[ Cell::iS_tmp ];
		Solution       &S_out         = cell.solutions[ Cell::iS_out ];
		cell.steadyState();
		
		
		
		
		//======================================================================
		//
		// compute time constants
		//
		//======================================================================
		const Real omega_AE  = cell.pAE->omega( S_in["HCO3-"], S_in["H+"]);
		const Real omega_NHE = cell.pNHE->omega( S_in["H+"] );
		std::cerr << "-- omega_AE="  << omega_AE  << std::endl;
		std::cerr << "-- omega_NHE=" << omega_NHE << std::endl;
		const Real omega_relax = -(omega_AE-omega_NHE);
		std::cerr << "-- relax at " << omega_relax << "/s" << std::endl;
		
		const Real      zeta0    = cell.Em * cell.F_over_RT;
		const Real      omega_K  = cell.library["K+"] .omega(zeta0,cell.Volume);
		const Real      omega_Na = cell.library["Na+"].omega(zeta0,cell.Volume);
		const Real      omega_Cl = cell.library["Cl-"].omega(zeta0,cell.Volume);
		
		std::cerr << "-- zeta0    = " << zeta0 << std::endl;
		std::cerr << "-- omega_K  = " << omega_K  << " /s"    << std::endl;
		std::cerr << "-- omega_Na = " << omega_Na << " /s"    << std::endl;
		std::cerr << "-- omega_Cl = " << omega_Cl << " /s"    << std::endl;
		const Real sum_dLeak_dzeta = S_tmp.set_dLeak_dzeta(zeta0, cell.Volume, S_in,S_out);
		std::cerr << "-- dLeak/dzeta=" << sum_dLeak_dzeta << std::endl;
		const Real drhoNaK_dzeta   = cell.pNaK->drho_dzeta(zeta0, cell);
		std::cerr << "-- drhoNak/dzeta=" << drhoNaK_dzeta << std::endl;
		const Real dincrease_dzeta = sum_dLeak_dzeta - drhoNaK_dzeta;
		//std::cerr << "-- dincrease/dzeta=" << dincrease_dzeta << std::endl;
		const Real Em_fac   =  (Faraday * cell.Volume * 1e-15)/cell.Capa;
		const Real omega_Em = cell.F_over_RT * Em_fac * dincrease_dzeta; //! volume: microns^3 => Litre;
		std::cerr << "-- omega_Em=" << omega_Em << " /s" << std::endl;
		
		matrix<Real> Omega(5,5);
		
		//-- variatons of Em
		/* /Em */ Omega[1][1] =  omega_Em;
		/* /K  */ Omega[1][2] =  Em_fac * omega_K;
		/* /Na */ Omega[1][3] =  Em_fac * omega_Na;
		/* /Cl */ Omega[1][4] = -Em_fac * omega_Cl;
		
		//-- variations of K
		/* /Em */ Omega[2][1] = cell.F_over_RT * (S_tmp[ "K+" ]+2*drhoNaK_dzeta);
		/* /K  */ Omega[2][2] = omega_K;
		
		//-- variations of Na
		/* /Em */ Omega[3][1] = cell.F_over_RT * (S_tmp[ "Na+" ]-3*drhoNaK_dzeta);
		/* /Na */ Omega[3][3] = omega_Na;
		/* /H+ */ Omega[3][5] = omega_NHE;
		
		//-- variations of Cl
		/* /Em */ Omega[4][1] = cell.F_over_RT * (S_tmp[ "Cl-" ] );
		/* /Cl */ Omega[4][4] = omega_Cl;
		/* /H  */ Omega[4][5] = omega_AE;
		
		//-- variations of H
		/* /H */ Omega[5][5] = omega_AE-omega_NHE; 
		
		std::cerr << "Omega=" << Omega << std::endl;
		
		std::cerr.unsetf( std::ios::floatfield );
		std::cerr << "tau_Em = "   << -1.0 / Omega[1][1]        << " s" << std::endl;
		std::cerr << "tau_K  = "   << (-1.0 / Omega[2][2])/60.0 << " mn" << std::endl;
		std::cerr << "tau_Na = "   << (-1.0 / Omega[3][3])/60.0 << " mn" << std::endl;
		std::cerr << "tau_Cl = "   << (-1.0 / Omega[4][4])/60.0 << " mn" << std::endl;

		
		//======================================================================
		//
		// Special values for H
		//
		//======================================================================
		//cell.simpleProtonRate(L);
		
		//======================================================================
		//
		// prepare simulation
		//
		//======================================================================
		Real   t    = 0;
		size_t iter = 0;
		Driver::Type odeint(1e-4);
		Driver::Init( odeint, cell );
        
        std::cerr << "Intrinsic DeltaQ:  " << cell.DeltaQ() << std::endl;
        std::cerr << "Intrinsic DeltaEm: " << cell.DeltaEm() << std::endl;
        std::cerr << "Initial   Charge:  " << cell.Em << std::endl;
        std::cerr << "Excess of Charge:  " << cell.excessConc( cell.Em ) << " M" << std::endl;
        
        
        std::cerr << std::endl;
		std::cerr << "Running: dt=" << dt << ", tmax=" << tmax << std::endl;
		
		
		Lua::Function<Real> P_CO2(L,"P_CO2");
		//----------------------------------------------------------------------
		// save starting point, with header
		//----------------------------------------------------------------------
		{
			ios::ocstream fp(output_name,false);
			fp("#t(mn) Em(mV) pH");
			S_in.write_header(fp);
			fp(" PCO2(mmHg)");
			fp("\n");
			save_data(fp,t,cell,P_CO2(t));
		}
		
#if defined(SAVE_LEAK)
		const string leak_name   = "leak-" + output_name;
		{
			ios::ocstream fp(leak_name,false);
			fp("#t(mn) Em(mV) pH");
			S_in.write_header(fp);
			fp("\n");
			save_leak(fp, t, cell);
		}
#endif
		const Real Q0 = S_in.sumCharge();
		std::cerr.flush();
		for(;;)
		{
			++iter;
			const Real t_old = (iter-1) * dt;
			t                = (iter)   * dt;
			Driver::Step( odeint, cell, t_old, t );		
			const Real dQ= S_in.sumCharge() - Q0;
			//------------------------------------------------------------------
			// save data
			//------------------------------------------------------------------
			fprintf( stderr, "t= %6.1f   , P_CO2=%10.6f, dQ=%7.2e   \r", t, P_CO2(t)*760.0, dQ);
			fflush( stderr );
			{
				ios::ocstream fp(output_name,true);
				save_data(fp,t,cell,P_CO2(t));
			}
#if defined(SAVE_LEAK)
			{
				ios::ocstream fp(leak_name,true);
				save_leak(fp, t, cell);
				
			}
#endif
			if( t >= tmax )
				break;
		}
		fprintf( stderr, "\n");
		fflush( stderr );
		
		return 0;
	}
	catch( const exception &e )
	{
		std::cerr << std::endl;
		std::cerr << e.when() << std::endl;
		std::cerr << e.what() << std::endl;
	}
	catch(...)
	{
		std::cerr << std::endl;
		std::cerr << "Unhandled Exception in " << progname << std::endl;
	}
	return -1;
}
